############
############
############	This file creates the true population of 1,000,000 people living in
############	25 different units. The first ten are large units (7%) and the latter
############	15 units are small (3%).
############	The individual variables (normally income, age, education, ...) is modeled 
############	here as gamma1, gamma2, gamma3.
############
############	Lucas Leemann, lleemann@gmail.com
############	December 2015 
############
############
############
############
############









############
############ Be aware that this simulation can take very long - it took 18 hours on a MacBook Pro (2.8 GHz, 16 GB) 
############










library(MASS)
library(lme4)


# set the working directory
setwd(".../Final replication file")



# Generate Population Data
source("Draws rho 0 vers 1.R")
source("Draws rho 0 vers B04.R")
source("Draws rho 0 vers B085.R")
source("Draws rho 0.2 vers 1.R")
source("Draws rho 0.2 vers B04.R")
source("Draws rho 0.2 vers B085.R")
source("Draws rho 0.4 vers 1.R")
source("Draws rho 0.4 vers B04.R")
source("Draws rho 0.4 vers B085.R")
source("Draws rho 0.6 vers 1.R")
source("Draws rho 0.6 vers B04.R")
source("Draws rho 0.6 vers B085.R")


# Run MC Analysis

source("MC A 0500 0.0 B00.R")
source("MC A 0500 0.0 B04.R")
source("MC A 0500 0.0 B085.R")
source("MC A 0500 0.2 B00.R")
source("MC A 0500 0.2 B04.R")
source("MC A 0500 0.2 B085.R")
source("MC A 0500 0.4 B00.R")
source("MC A 0500 0.4 B04.R")
source("MC A 0500 0.4 B085.R")
source("MC A 0500 0.6 B00.R")
source("MC A 0500 0.6 B04.R")
source("MC A 0500 0.6 B085.R")


source("MC A 1000 0.0 B00.R")
source("MC A 1000 0.0 B04.R")
source("MC A 1000 0.0 B085.R")
source("MC A 1000 0.2 B00.R")
source("MC A 1000 0.2 B04.R")
source("MC A 1000 0.2 B085.R")
source("MC A 1000 0.4 B00.R")
source("MC A 1000 0.4 B04.R")
source("MC A 1000 0.4 B085.R")
source("MC A 1000 0.6 B00.R")
source("MC A 1000 0.6 B04.R")
source("MC A 1000 0.6 B085.R")

source("MC A 2000 0.0 B00.R")
source("MC A 2000 0.0 B04.R")
source("MC A 2000 0.0 B085.R")
source("MC A 2000 0.2 B00.R")
source("MC A 2000 0.2 B04.R")
source("MC A 2000 0.2 B085.R")
source("MC A 2000 0.4 B00.R")
source("MC A 2000 0.4 B04.R")
source("MC A 2000 0.4 B085.R")
source("MC A 2000 0.6 B00.R")
source("MC A 2000 0.6 B04.R")
source("MC A 2000 0.6 B085.R")




### Illustrate MC Simulations



load("sim.0500.00.B00.RData")
load("sim.1000.00.B00.RData")
load("sim.2000.00.B00.RData")
load("sim.0500.02.B00.RData")
load("sim.1000.02.B00.RData")
load("sim.2000.02.B00.RData")
load("sim.0500.04.B00.RData")
load("sim.1000.04.B00.RData")
load("sim.2000.04.B00.RData")
load("sim.0500.06.B00.RData")
load("sim.1000.06.B00.RData")
load("sim.2000.06.B00.RData")

load("sim.0500.00.B04.RData")
load("sim.1000.00.B04.RData")
load("sim.2000.00.B04.RData")
load("sim.0500.02.B04.RData")
load("sim.1000.02.B04.RData")
load("sim.2000.02.B04.RData")
load("sim.0500.04.B04.RData")
load("sim.1000.04.B04.RData")
load("sim.2000.04.B04.RData")
load("sim.0500.06.B04.RData")
load("sim.1000.06.B04.RData")
load("sim.2000.06.B04.RData")

load("sim.0500.00.B085.RData")
load("sim.1000.00.B085.RData")
load("sim.2000.00.B085.RData")
load("sim.0500.02.B085.RData")
load("sim.1000.02.B085.RData")
load("sim.2000.02.B085.RData")
load("sim.0500.04.B085.RData")
load("sim.1000.04.B085.RData")
load("sim.2000.04.B085.RData")
load("sim.0500.06.B085.RData")
load("sim.1000.06.B085.RData")
load("sim.2000.06.B085.RData")

matrix.mse.B00 <- array(NA,c(100,12,3))

matrix.mse.B00[,1,1] <- rowMeans(abs(sim.0500.0.0rho[8][[1]]))			# mse Joint
matrix.mse.B00[,1,2] <- rowMeans(abs(sim.0500.0.0rho[9][[1]]))			# mse Marginal
matrix.mse.B00[,1,3] <- rowMeans(abs(sim.0500.0.0rho[10][[1]]))			# mse Disaggregation

matrix.mse.B00[,2,1] <- rowMeans(abs(sim.0500.0.2rho[8][[1]]))			# mse Joint
matrix.mse.B00[,2,2] <- rowMeans(abs(sim.0500.0.2rho[9][[1]]))			# mse Marginal
matrix.mse.B00[,2,3] <- rowMeans(abs(sim.0500.0.2rho[10][[1]]))			# mse Disaggregation

matrix.mse.B00[,3,1] <- rowMeans(abs(sim.0500.0.4rho[8][[1]]))			# mse Joint
matrix.mse.B00[,3,2] <- rowMeans(abs(sim.0500.0.4rho[9][[1]]))			# mse Marginal
matrix.mse.B00[,3,3] <- rowMeans(abs(sim.0500.0.4rho[10][[1]]))			# mse Disaggregation

matrix.mse.B00[,4,1] <- rowMeans(abs(sim.0500.0.6rho[8][[1]]))			# mse Joint
matrix.mse.B00[,4,2] <- rowMeans(abs(sim.0500.0.6rho[9][[1]]))			# mse Marginal
matrix.mse.B00[,4,3] <- rowMeans(abs(sim.0500.0.6rho[10][[1]]))			# mse Disaggregation

matrix.mse.B00[,5,1] <- rowMeans(abs(sim.1000.0.0rho[8][[1]]))			# mse Joint
matrix.mse.B00[,5,2] <- rowMeans(abs(sim.1000.0.0rho[9][[1]]))			# mse Marginal
matrix.mse.B00[,5,3] <- rowMeans(abs(sim.1000.0.0rho[10][[1]]))			# mse Disaggregation

matrix.mse.B00[,6,1] <- rowMeans(abs(sim.1000.0.2rho[8][[1]]))			# mse Joint
matrix.mse.B00[,6,2] <- rowMeans(abs(sim.1000.0.2rho[9][[1]]))			# mse Marginal
matrix.mse.B00[,6,3] <- rowMeans(abs(sim.1000.0.2rho[10][[1]]))			# mse Disaggregation

matrix.mse.B00[,7,1] <- rowMeans(abs(sim.1000.0.4rho[8][[1]]))			# mse Joint
matrix.mse.B00[,7,2] <- rowMeans(abs(sim.1000.0.4rho[9][[1]]))			# mse Marginal
matrix.mse.B00[,7,3] <- rowMeans(abs(sim.1000.0.4rho[10][[1]]))			# mse Disaggregation

matrix.mse.B00[,8,1] <- rowMeans(abs(sim.1000.0.6rho[8][[1]]))			# mse Joint
matrix.mse.B00[,8,2] <- rowMeans(abs(sim.1000.0.6rho[9][[1]]))			# mse Marginal
matrix.mse.B00[,8,3] <- rowMeans(abs(sim.1000.0.6rho[10][[1]]))			# mse Disaggregation

matrix.mse.B00[,9,1] <- rowMeans(abs(sim.2000.0.0rho[8][[1]]))			# mse Joint
matrix.mse.B00[,9,2] <- rowMeans(abs(sim.2000.0.0rho[9][[1]]))			# mse Marginal
matrix.mse.B00[,9,3] <- rowMeans(abs(sim.2000.0.0rho[10][[1]]))			# mse Disaggregation

matrix.mse.B00[,10,1] <- rowMeans(abs(sim.2000.0.2rho[8][[1]]))			# mse Joint
matrix.mse.B00[,10,2] <- rowMeans(abs(sim.2000.0.2rho[9][[1]]))			# mse Marginal
matrix.mse.B00[,10,3] <- rowMeans(abs(sim.2000.0.2rho[10][[1]]))			# mse Disaggregation

matrix.mse.B00[,11,1] <- rowMeans(abs(sim.2000.0.4rho[8][[1]]))			# mse Joint
matrix.mse.B00[,11,2] <- rowMeans(abs(sim.2000.0.4rho[9][[1]]))			# mse Marginal
matrix.mse.B00[,11,3] <- rowMeans(abs(sim.2000.0.4rho[10][[1]]))			# mse Disaggregation

matrix.mse.B00[,12,1] <- rowMeans(abs(sim.2000.0.6rho[8][[1]]))		# mse Joint
matrix.mse.B00[,12,2] <- rowMeans(abs(sim.2000.0.6rho[9][[1]]))		# mse Marginal
matrix.mse.B00[,12,3] <- rowMeans(abs(sim.2000.0.6rho[10][[1]]))		# mse Disaggregation


matrix.mse.B04 <- array(NA,c(100,12,3))

matrix.mse.B04[,1,1] <- rowMeans(abs(sim.0500.0.0rhoB04[8][[1]]))			# mse Joint
matrix.mse.B04[,1,2] <- rowMeans(abs(sim.0500.0.0rhoB04[9][[1]]))			# mse Marginal
matrix.mse.B04[,1,3] <- rowMeans(abs(sim.0500.0.0rhoB04[10][[1]]))			# mse Disaggregation

matrix.mse.B04[,2,1] <- rowMeans(abs(sim.0500.0.2rhoB04[8][[1]]))			# mse Joint
matrix.mse.B04[,2,2] <- rowMeans(abs(sim.0500.0.2rhoB04[9][[1]]))			# mse Marginal
matrix.mse.B04[,2,3] <- rowMeans(abs(sim.0500.0.2rhoB04[10][[1]]))			# mse Disaggregation

matrix.mse.B04[,3,1] <- rowMeans(abs(sim.0500.0.4rhoB04[8][[1]]))			# mse Joint
matrix.mse.B04[,3,2] <- rowMeans(abs(sim.0500.0.4rhoB04[9][[1]]))			# mse Marginal
matrix.mse.B04[,3,3] <- rowMeans(abs(sim.0500.0.4rhoB04[10][[1]]))			# mse Disaggregation

matrix.mse.B04[,4,1] <- rowMeans(abs(sim.0500.0.6rhoB04[8][[1]]))			# mse Joint
matrix.mse.B04[,4,2] <- rowMeans(abs(sim.0500.0.6rhoB04[9][[1]]))			# mse Marginal
matrix.mse.B04[,4,3] <- rowMeans(abs(sim.0500.0.6rhoB04[10][[1]]))			# mse Disaggregation

matrix.mse.B04[,5,1] <- rowMeans(abs(sim.1000.0.0rhoB04[8][[1]]))			# mse Joint
matrix.mse.B04[,5,2] <- rowMeans(abs(sim.1000.0.0rhoB04[9][[1]]))			# mse Marginal
matrix.mse.B04[,5,3] <- rowMeans(abs(sim.1000.0.0rhoB04[10][[1]]))			# mse Disaggregation

matrix.mse.B04[,6,1] <- rowMeans(abs(sim.1000.0.2rhoB04[8][[1]]))			# mse Joint
matrix.mse.B04[,6,2] <- rowMeans(abs(sim.1000.0.2rhoB04[9][[1]]))			# mse Marginal
matrix.mse.B04[,6,3] <- rowMeans(abs(sim.1000.0.2rhoB04[10][[1]]))			# mse Disaggregation

matrix.mse.B04[,7,1] <- rowMeans(abs(sim.1000.0.4rhoB04[8][[1]]))			# mse Joint
matrix.mse.B04[,7,2] <- rowMeans(abs(sim.1000.0.4rhoB04[9][[1]]))			# mse Marginal
matrix.mse.B04[,7,3] <- rowMeans(abs(sim.1000.0.4rhoB04[10][[1]]))			# mse Disaggregation

matrix.mse.B04[,8,1] <- rowMeans(abs(sim.1000.0.6rhoB04[8][[1]]))			# mse Joint
matrix.mse.B04[,8,2] <- rowMeans(abs(sim.1000.0.6rhoB04[9][[1]]))			# mse Marginal
matrix.mse.B04[,8,3] <- rowMeans(abs(sim.1000.0.6rhoB04[10][[1]]))			# mse Disaggregation

matrix.mse.B04[,9,1] <- rowMeans(abs(sim.2000.0.0rhoB04[8][[1]]))			# mse Joint
matrix.mse.B04[,9,2] <- rowMeans(abs(sim.2000.0.0rhoB04[9][[1]]))			# mse Marginal
matrix.mse.B04[,9,3] <- rowMeans(abs(sim.2000.0.0rhoB04[10][[1]]))			# mse Disaggregation

matrix.mse.B04[,10,1] <- rowMeans(abs(sim.2000.0.2rhoB04[8][[1]]))			# mse Joint
matrix.mse.B04[,10,2] <- rowMeans(abs(sim.2000.0.2rhoB04[9][[1]]))			# mse Marginal
matrix.mse.B04[,10,3] <- rowMeans(abs(sim.2000.0.2rhoB04[10][[1]]))			# mse Disaggregation

matrix.mse.B04[,11,1] <- rowMeans(abs(sim.2000.0.4rhoB04[8][[1]]))			# mse Joint
matrix.mse.B04[,11,2] <- rowMeans(abs(sim.2000.0.4rhoB04[9][[1]]))			# mse Marginal
matrix.mse.B04[,11,3] <- rowMeans(abs(sim.2000.0.4rhoB04[10][[1]]))			# mse Disaggregation

matrix.mse.B04[,12,1] <- rowMeans(abs(sim.2000.0.6rhoB04[8][[1]]))		# mse Joint
matrix.mse.B04[,12,2] <- rowMeans(abs(sim.2000.0.6rhoB04[9][[1]]))		# mse Marginal
matrix.mse.B04[,12,3] <- rowMeans(abs(sim.2000.0.6rhoB04[10][[1]]))		# mse Disaggregation

matrix.mse.B085 <- array(NA,c(100,12,3))

matrix.mse.B085[,1,1] <- rowMeans(abs(sim.0500.0.0rhoB085[8][[1]]))			# mse Joint
matrix.mse.B085[,1,2] <- rowMeans(abs(sim.0500.0.0rhoB085[9][[1]]))			# mse Marginal
matrix.mse.B085[,1,3] <- rowMeans(abs(sim.0500.0.0rhoB085[10][[1]]))			# mse Disaggregation

matrix.mse.B085[,2,1] <- rowMeans(abs(sim.0500.0.2rhoB085[8][[1]]))			# mse Joint
matrix.mse.B085[,2,2] <- rowMeans(abs(sim.0500.0.2rhoB085[9][[1]]))			# mse Marginal
matrix.mse.B085[,2,3] <- rowMeans(abs(sim.0500.0.2rhoB085[10][[1]]))			# mse Disaggregation

matrix.mse.B085[,3,1] <- rowMeans(abs(sim.0500.0.4rhoB085[8][[1]]))			# mse Joint
matrix.mse.B085[,3,2] <- rowMeans(abs(sim.0500.0.4rhoB085[9][[1]]))			# mse Marginal
matrix.mse.B085[,3,3] <- rowMeans(abs(sim.0500.0.4rhoB085[10][[1]]))			# mse Disaggregation

matrix.mse.B085[,4,1] <- rowMeans(abs(sim.0500.0.6rhoB085[8][[1]]))			# mse Joint
matrix.mse.B085[,4,2] <- rowMeans(abs(sim.0500.0.6rhoB085[9][[1]]))			# mse Marginal
matrix.mse.B085[,4,3] <- rowMeans(abs(sim.0500.0.6rhoB085[10][[1]]))			# mse Disaggregation

matrix.mse.B085[,5,1] <- rowMeans(abs(sim.1000.0.0rhoB085[8][[1]]))			# mse Joint
matrix.mse.B085[,5,2] <- rowMeans(abs(sim.1000.0.0rhoB085[9][[1]]))			# mse Marginal
matrix.mse.B085[,5,3] <- rowMeans(abs(sim.1000.0.0rhoB085[10][[1]]))			# mse Disaggregation

matrix.mse.B085[,6,1] <- rowMeans(abs(sim.1000.0.2rhoB085[8][[1]]))			# mse Joint
matrix.mse.B085[,6,2] <- rowMeans(abs(sim.1000.0.2rhoB085[9][[1]]))			# mse Marginal
matrix.mse.B085[,6,3] <- rowMeans(abs(sim.1000.0.2rhoB085[10][[1]]))			# mse Disaggregation

matrix.mse.B085[,7,1] <- rowMeans(abs(sim.1000.0.4rhoB085[8][[1]]))			# mse Joint
matrix.mse.B085[,7,2] <- rowMeans(abs(sim.1000.0.4rhoB085[9][[1]]))			# mse Marginal
matrix.mse.B085[,7,3] <- rowMeans(abs(sim.1000.0.4rhoB085[10][[1]]))			# mse Disaggregation

matrix.mse.B085[,8,1] <- rowMeans(abs(sim.1000.0.6rhoB085[8][[1]]))			# mse Joint
matrix.mse.B085[,8,2] <- rowMeans(abs(sim.1000.0.6rhoB085[9][[1]]))			# mse Marginal
matrix.mse.B085[,8,3] <- rowMeans(abs(sim.1000.0.6rhoB085[10][[1]]))			# mse Disaggregation

matrix.mse.B085[,9,1] <- rowMeans(abs(sim.2000.0.0rhoB085[8][[1]]))			# mse Joint
matrix.mse.B085[,9,2] <- rowMeans(abs(sim.2000.0.0rhoB085[9][[1]]))			# mse Marginal
matrix.mse.B085[,9,3] <- rowMeans(abs(sim.2000.0.0rhoB085[10][[1]]))			# mse Disaggregation

matrix.mse.B085[,10,1] <- rowMeans(abs(sim.2000.0.2rhoB085[8][[1]]))			# mse Joint
matrix.mse.B085[,10,2] <- rowMeans(abs(sim.2000.0.2rhoB085[9][[1]]))			# mse Marginal
matrix.mse.B085[,10,3] <- rowMeans(abs(sim.2000.0.2rhoB085[10][[1]]))			# mse Disaggregation

matrix.mse.B085[,11,1] <- rowMeans(abs(sim.2000.0.4rhoB085[8][[1]]))			# mse Joint
matrix.mse.B085[,11,2] <- rowMeans(abs(sim.2000.0.4rhoB085[9][[1]]))			# mse Marginal
matrix.mse.B085[,11,3] <- rowMeans(abs(sim.2000.0.4rhoB085[10][[1]]))			# mse Disaggregation

matrix.mse.B085[,12,1] <- rowMeans(abs(sim.2000.0.6rhoB085[8][[1]]))		# mse Joint
matrix.mse.B085[,12,2] <- rowMeans(abs(sim.2000.0.6rhoB085[9][[1]]))		# mse Marginal
matrix.mse.B085[,12,3] <- rowMeans(abs(sim.2000.0.6rhoB085[10][[1]]))		# mse Disaggregation


### Plot Simple

LL.dot <- function(x, quant){
		x <- sort(x)	
		N <- length(x)
		low <- floor(quant[1]*N)
		high <- ceiling(quant[2]*N)
		outL <- rep(NA, 3)
		outL[1] <- x[low]
		outL[2] <- mean(x)
		outL[3] <- x[high]
		print(outL)
	}

# example
LL.dot(matrix.mse.B00[,1,1],quant=c(0.025,0.975))
qq <- c(0.025,0.975)

#prepare objects
# sample: 500
mae.jo.0500.B00 <- LL.dot(matrix.mse.B00[,1,1], quant=qq)
mae.ma.0500.B00 <- LL.dot(matrix.mse.B00[,1,2], quant=qq)
mae.di.0500.B00 <- LL.dot(matrix.mse.B00[,1,3], quant=qq)
mae.jo.0502.B00 <- LL.dot(matrix.mse.B00[,2,1], quant=qq)
mae.ma.0502.B00 <- LL.dot(matrix.mse.B00[,2,2], quant=qq)
mae.di.0502.B00 <- LL.dot(matrix.mse.B00[,2,3], quant=qq)
mae.jo.0504.B00 <- LL.dot(matrix.mse.B00[,3,1], quant=qq)
mae.ma.0504.B00 <- LL.dot(matrix.mse.B00[,3,2], quant=qq)
mae.di.0504.B00 <- LL.dot(matrix.mse.B00[,3,3], quant=qq)
mae.jo.0506.B00 <- LL.dot(matrix.mse.B00[,4,1], quant=qq)
mae.ma.0506.B00 <- LL.dot(matrix.mse.B00[,4,2], quant=qq)
mae.di.0506.B00 <- LL.dot(matrix.mse.B00[,4,3], quant=qq)
# sample: 1000
mae.jo.1000.B00 <- LL.dot(matrix.mse.B00[,5,1], quant=qq)
mae.ma.1000.B00 <- LL.dot(matrix.mse.B00[,5,2], quant=qq)
mae.di.1000.B00 <- LL.dot(matrix.mse.B00[,5,3], quant=qq)
mae.jo.1002.B00 <- LL.dot(matrix.mse.B00[,6,1], quant=qq)
mae.ma.1002.B00 <- LL.dot(matrix.mse.B00[,6,2], quant=qq)
mae.di.1002.B00 <- LL.dot(matrix.mse.B00[,6,3], quant=qq)
mae.jo.1004.B00 <- LL.dot(matrix.mse.B00[,7,1], quant=qq)
mae.ma.1004.B00 <- LL.dot(matrix.mse.B00[,7,2], quant=qq)
mae.di.1004.B00 <- LL.dot(matrix.mse.B00[,7,3], quant=qq)
mae.jo.1006.B00 <- LL.dot(matrix.mse.B00[,8,1], quant=qq)
mae.ma.1006.B00 <- LL.dot(matrix.mse.B00[,8,2], quant=qq)
mae.di.1006.B00 <- LL.dot(matrix.mse.B00[,8,3], quant=qq)
# sample: 2000
mae.jo.2000.B00 <- LL.dot(matrix.mse.B00[,9,1], quant=qq)
mae.ma.2000.B00 <- LL.dot(matrix.mse.B00[,9,2], quant=qq)
mae.di.2000.B00 <- LL.dot(matrix.mse.B00[,9,3], quant=qq)
mae.jo.2002.B00 <- LL.dot(matrix.mse.B00[,10,1], quant=qq)
mae.ma.2002.B00 <- LL.dot(matrix.mse.B00[,10,2], quant=qq)
mae.di.2002.B00 <- LL.dot(matrix.mse.B00[,10,3], quant=qq)
mae.jo.2004.B00 <- LL.dot(matrix.mse.B00[,11,1], quant=qq)
mae.ma.2004.B00 <- LL.dot(matrix.mse.B00[,11,2], quant=qq)
mae.di.2004.B00 <- LL.dot(matrix.mse.B00[,11,3], quant=qq)
mae.jo.2006.B00 <- LL.dot(matrix.mse.B00[,12,1], quant=qq)
mae.ma.2006.B00 <- LL.dot(matrix.mse.B00[,12,2], quant=qq)
mae.di.2006.B00 <- LL.dot(matrix.mse.B00[,12,3], quant=qq)

# sample: 500
mae.jo.0500.B04 <- LL.dot(matrix.mse.B04[,1,1], quant=qq)
mae.ma.0500.B04 <- LL.dot(matrix.mse.B04[,1,2], quant=qq)
mae.di.0500.B04 <- LL.dot(matrix.mse.B04[,1,3], quant=qq)
mae.jo.0502.B04 <- LL.dot(matrix.mse.B04[,2,1], quant=qq)
mae.ma.0502.B04 <- LL.dot(matrix.mse.B04[,2,2], quant=qq)
mae.di.0502.B04 <- LL.dot(matrix.mse.B04[,2,3], quant=qq)
mae.jo.0504.B04 <- LL.dot(matrix.mse.B04[,3,1], quant=qq)
mae.ma.0504.B04 <- LL.dot(matrix.mse.B04[,3,2], quant=qq)
mae.di.0504.B04 <- LL.dot(matrix.mse.B04[,3,3], quant=qq)
mae.jo.0506.B04 <- LL.dot(matrix.mse.B04[,4,1], quant=qq)
mae.ma.0506.B04 <- LL.dot(matrix.mse.B04[,4,2], quant=qq)
mae.di.0506.B04 <- LL.dot(matrix.mse.B04[,4,3], quant=qq)
# sample: 1000
mae.jo.1000.B04 <- LL.dot(matrix.mse.B04[,5,1], quant=qq)
mae.ma.1000.B04 <- LL.dot(matrix.mse.B04[,5,2], quant=qq)
mae.di.1000.B04 <- LL.dot(matrix.mse.B04[,5,3], quant=qq)
mae.jo.1002.B04 <- LL.dot(matrix.mse.B04[,6,1], quant=qq)
mae.ma.1002.B04 <- LL.dot(matrix.mse.B04[,6,2], quant=qq)
mae.di.1002.B04 <- LL.dot(matrix.mse.B04[,6,3], quant=qq)
mae.jo.1004.B04 <- LL.dot(matrix.mse.B04[,7,1], quant=qq)
mae.ma.1004.B04 <- LL.dot(matrix.mse.B04[,7,2], quant=qq)
mae.di.1004.B04 <- LL.dot(matrix.mse.B04[,7,3], quant=qq)
mae.jo.1006.B04 <- LL.dot(matrix.mse.B04[,8,1], quant=qq)
mae.ma.1006.B04 <- LL.dot(matrix.mse.B04[,8,2], quant=qq)
mae.di.1006.B04 <- LL.dot(matrix.mse.B04[,8,3], quant=qq)
# sample: 2000
mae.jo.2000.B04 <- LL.dot(matrix.mse.B04[,9,1], quant=qq)
mae.ma.2000.B04 <- LL.dot(matrix.mse.B04[,9,2], quant=qq)
mae.di.2000.B04 <- LL.dot(matrix.mse.B04[,9,3], quant=qq)
mae.jo.2002.B04 <- LL.dot(matrix.mse.B04[,10,1], quant=qq)
mae.ma.2002.B04 <- LL.dot(matrix.mse.B04[,10,2], quant=qq)
mae.di.2002.B04 <- LL.dot(matrix.mse.B04[,10,3], quant=qq)
mae.jo.2004.B04 <- LL.dot(matrix.mse.B04[,11,1], quant=qq)
mae.ma.2004.B04 <- LL.dot(matrix.mse.B04[,11,2], quant=qq)
mae.di.2004.B04 <- LL.dot(matrix.mse.B04[,11,3], quant=qq)
mae.jo.2006.B04 <- LL.dot(matrix.mse.B04[,12,1], quant=qq)
mae.ma.2006.B04 <- LL.dot(matrix.mse.B04[,12,2], quant=qq)
mae.di.2006.B04 <- LL.dot(matrix.mse.B04[,12,3], quant=qq)


# sample: 500
mae.jo.0500.B085 <- LL.dot(matrix.mse.B085[,1,1], quant=qq)
mae.ma.0500.B085 <- LL.dot(matrix.mse.B085[,1,2], quant=qq)
mae.di.0500.B085 <- LL.dot(matrix.mse.B085[,1,3], quant=qq)
mae.jo.0502.B085 <- LL.dot(matrix.mse.B085[,2,1], quant=qq)
mae.ma.0502.B085 <- LL.dot(matrix.mse.B085[,2,2], quant=qq)
mae.di.0502.B085 <- LL.dot(matrix.mse.B085[,2,3], quant=qq)
mae.jo.0504.B085 <- LL.dot(matrix.mse.B085[,3,1], quant=qq)
mae.ma.0504.B085 <- LL.dot(matrix.mse.B085[,3,2], quant=qq)
mae.di.0504.B085 <- LL.dot(matrix.mse.B085[,3,3], quant=qq)
mae.jo.0506.B085 <- LL.dot(matrix.mse.B085[,4,1], quant=qq)
mae.ma.0506.B085 <- LL.dot(matrix.mse.B085[,4,2], quant=qq)
mae.di.0506.B085 <- LL.dot(matrix.mse.B085[,4,3], quant=qq)
# sample: 1000
mae.jo.1000.B085 <- LL.dot(matrix.mse.B085[,5,1], quant=qq)
mae.ma.1000.B085 <- LL.dot(matrix.mse.B085[,5,2], quant=qq)
mae.di.1000.B085 <- LL.dot(matrix.mse.B085[,5,3], quant=qq)
mae.jo.1002.B085 <- LL.dot(matrix.mse.B085[,6,1], quant=qq)
mae.ma.1002.B085 <- LL.dot(matrix.mse.B085[,6,2], quant=qq)
mae.di.1002.B085 <- LL.dot(matrix.mse.B085[,6,3], quant=qq)
mae.jo.1004.B085 <- LL.dot(matrix.mse.B085[,7,1], quant=qq)
mae.ma.1004.B085 <- LL.dot(matrix.mse.B085[,7,2], quant=qq)
mae.di.1004.B085 <- LL.dot(matrix.mse.B085[,7,3], quant=qq)
mae.jo.1006.B085 <- LL.dot(matrix.mse.B085[,8,1], quant=qq)
mae.ma.1006.B085 <- LL.dot(matrix.mse.B085[,8,2], quant=qq)
mae.di.1006.B085 <- LL.dot(matrix.mse.B085[,8,3], quant=qq)
# sample: 2000
mae.jo.2000.B085 <- LL.dot(matrix.mse.B085[,9,1], quant=qq)
mae.ma.2000.B085 <- LL.dot(matrix.mse.B085[,9,2], quant=qq)
mae.di.2000.B085 <- LL.dot(matrix.mse.B085[,9,3], quant=qq)
mae.jo.2002.B085 <- LL.dot(matrix.mse.B085[,10,1], quant=qq)
mae.ma.2002.B085 <- LL.dot(matrix.mse.B085[,10,2], quant=qq)
mae.di.2002.B085 <- LL.dot(matrix.mse.B085[,10,3], quant=qq)
mae.jo.2004.B085 <- LL.dot(matrix.mse.B085[,11,1], quant=qq)
mae.ma.2004.B085 <- LL.dot(matrix.mse.B085[,11,2], quant=qq)
mae.di.2004.B085 <- LL.dot(matrix.mse.B085[,11,3], quant=qq)
mae.jo.2006.B085 <- LL.dot(matrix.mse.B085[,12,1], quant=qq)
mae.ma.2006.B085 <- LL.dot(matrix.mse.B085[,12,2], quant=qq)
mae.di.2006.B085 <- LL.dot(matrix.mse.B085[,12,3], quant=qq)

### New Plot
png("Figure01.png", width = 1200, height = 1400, points=18)
par(mfrow=c(4,3),mar=c(4,6,3,1), oma=c(1,2,0,1))
#par(mfrow=c(4,3), family="CMU Serif",mar=c(4,6,3,1), oma=c(1,2,0,1))

segmentsL <- function(x1,x2,x3,x4, widL, col, lwd, lty){
	segments(x1,x2,x3,x4, col=col, lwd=lwd, lty=lty)
	segments(x1-widL/2,x2,x1+widL/2,x2,  col=col, lwd=lwd, lty=lty)
	segments(x1-widL/2,x4,x1+widL/2,x4,  col=col, lwd=lwd, lty=lty)
	}


q <- 0
plot(0.98, mae.jo.0500.B00[2], ylim=c(0,0.15), xlim=c(0.8,4.2), xaxt="n", cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3, ylab=expression(paste("Average support 50%, ", delta[c1], "=0.0")), cex.lab=1.5, cex.axis=1.5, cex.main=1.8, main="Sample Size: 500", xlab="")
	segmentsL(0.98,mae.jo.0500.B00[1],0.98,mae.jo.0500.B00[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02, mae.ma.0500.B00[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02,mae.ma.0500.B00[1],1.02,mae.ma.0500.B00[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1, mae.di.0500.B00[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1,mae.di.0500.B00[1],1,mae.di.0500.B00[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 1
	points(0.98+q, mae.jo.0502.B00[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.0502.B00[1],0.98+q,mae.jo.0502.B00[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.0502.B00[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.0502.B00[1],1.02+q,mae.ma.0502.B00[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.0502.B00[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.0502.B00[1],1+q,mae.di.0502.B00[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 2
	points(0.98+q, mae.jo.0504.B00[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.0504.B00[1],0.98+q,mae.jo.0504.B00[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.0504.B00[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.0504.B00[1],1.02+q,mae.ma.0504.B00[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.0504.B00[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.0504.B00[1],1+q,mae.di.0504.B00[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 3
	points(0.98+q, mae.jo.0506.B00[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.0506.B00[1],0.98+q,mae.jo.0506.B00[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.0506.B00[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.0506.B00[1],1.02+q,mae.ma.0506.B00[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.0506.B00[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.0506.B00[1],1+q,mae.di.0506.B00[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	# lines
	segments(1,mae.jo.0500.B00[2],2,mae.jo.0502.B00[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(2,mae.jo.0502.B00[2],3,mae.jo.0504.B00[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(3,mae.jo.0504.B00[2],4,mae.jo.0506.B00[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(1,mae.ma.0500.B00[2],2,mae.ma.0502.B00[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(2,mae.ma.0502.B00[2],3,mae.ma.0504.B00[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(3,mae.ma.0504.B00[2],4,mae.ma.0506.B00[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(1,mae.di.0500.B00[2],2,mae.di.0502.B00[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(2,mae.di.0502.B00[2],3,mae.di.0504.B00[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(3,mae.di.0504.B00[2],4,mae.di.0506.B00[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))




axis(1, at=c(1,2,3,4), labels=c(expression(paste(rho,"=0.0")), expression(paste(rho,"=0.2")), expression(paste(rho,"=0.4")), expression(paste(rho,"=0.6")) ), lwd.ticks=0, cex.axis=1.5)

# second plot
q <- 0
plot(0.98, mae.jo.1000.B00[2], ylim=c(0,0.15), xlim=c(0.8,4.2), xaxt="n", cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3, ylab="", cex.lab=1.8, cex.axis=1.5, cex.main=1.8, main="Sample Size: 1000", xlab="")
	segmentsL(0.98,mae.jo.1000.B00[1],0.98,mae.jo.1000.B00[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02, mae.ma.1000.B00[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02,mae.ma.1000.B00[1],1.02,mae.ma.1000.B00[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1, mae.di.1000.B00[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1,mae.di.1000.B00[1],1,mae.di.1000.B00[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 1
	points(0.98+q, mae.jo.1002.B00[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.1002.B00[1],0.98+q,mae.jo.1002.B00[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.1002.B00[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.1002.B00[1],1.02+q,mae.ma.1002.B00[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.1002.B00[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.1002.B00[1],1+q,mae.di.1002.B00[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 2
	points(0.98+q, mae.jo.1004.B00[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.1004.B00[1],0.98+q,mae.jo.1004.B00[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.1004.B00[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.1004.B00[1],1.02+q,mae.ma.1004.B00[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.1004.B00[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.1004.B00[1],1+q,mae.di.1004.B00[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 3
	points(0.98+q, mae.jo.1006.B00[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.1006.B00[1],0.98+q,mae.jo.1006.B00[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.1006.B00[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.1006.B00[1],1.02+q,mae.ma.1006.B00[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.1006.B00[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.1006.B00[1],1+q,mae.di.1006.B00[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	# lines
	segments(1,mae.jo.1000.B00[2],2,mae.jo.1002.B00[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(2,mae.jo.1002.B00[2],3,mae.jo.1004.B00[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(3,mae.jo.1004.B00[2],4,mae.jo.1006.B00[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(1,mae.ma.1000.B00[2],2,mae.ma.1002.B00[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(2,mae.ma.1002.B00[2],3,mae.ma.1004.B00[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(3,mae.ma.1004.B00[2],4,mae.ma.1006.B00[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(1,mae.di.1000.B00[2],2,mae.di.1002.B00[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(2,mae.di.1002.B00[2],3,mae.di.1004.B00[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(3,mae.di.1004.B00[2],4,mae.di.1006.B00[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))




axis(1, at=c(1,2,3,4), labels=c(expression(paste(rho,"=0.0")), expression(paste(rho,"=0.2")), expression(paste(rho,"=0.4")), expression(paste(rho,"=0.6")) ), lwd.ticks=0, cex.axis=1.5)



# third plot
q <- 0
plot(0.98, mae.jo.2000.B00[2], ylim=c(0,0.15), xlim=c(0.8,4.2), xaxt="n", cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3, ylab="", cex.lab=1.8, cex.axis=1.5, cex.main=1.8, main="Sample Size: 2000", xlab="")
	segmentsL(0.98,mae.jo.2000.B00[1],0.98,mae.jo.2000.B00[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02, mae.ma.2000.B00[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02,mae.ma.2000.B00[1],1.02,mae.ma.2000.B00[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1, mae.di.2000.B00[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1,mae.di.2000.B00[1],1,mae.di.2000.B00[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 1
	points(0.98+q, mae.jo.2002.B00[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.2002.B00[1],0.98+q,mae.jo.2002.B00[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.2002.B00[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.2002.B00[1],1.02+q,mae.ma.2002.B00[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.2002.B00[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.2002.B00[1],1+q,mae.di.2002.B00[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 2
	points(0.98+q, mae.jo.2004.B00[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.2004.B00[1],0.98+q,mae.jo.2004.B00[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.2004.B00[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.2004.B00[1],1.02+q,mae.ma.2004.B00[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.2004.B00[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.2004.B00[1],1+q,mae.di.2004.B00[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 3
	points(0.98+q, mae.jo.2006.B00[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.2006.B00[1],0.98+q,mae.jo.2006.B00[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.2006.B00[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.2006.B00[1],1.02+q,mae.ma.2006.B00[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.2006.B00[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.2006.B00[1],1+q,mae.di.2006.B00[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	# lines
	segments(1,mae.jo.2000.B00[2],2,mae.jo.2002.B00[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(2,mae.jo.2002.B00[2],3,mae.jo.2004.B00[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(3,mae.jo.2004.B00[2],4,mae.jo.2006.B00[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(1,mae.ma.2000.B00[2],2,mae.ma.2002.B00[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(2,mae.ma.2002.B00[2],3,mae.ma.2004.B00[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(3,mae.ma.2004.B00[2],4,mae.ma.2006.B00[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(1,mae.di.2000.B00[2],2,mae.di.2002.B00[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(2,mae.di.2002.B00[2],3,mae.di.2004.B00[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(3,mae.di.2004.B00[2],4,mae.di.2006.B00[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))




axis(1, at=c(1,2,3,4), labels=c(expression(paste(rho,"=0.0")), expression(paste(rho,"=0.2")), expression(paste(rho,"=0.4")), expression(paste(rho,"=0.6")) ), lwd.ticks=0, cex.axis=1.5)




q <- 0
plot(0.98, mae.jo.0500.B04[2], ylim=c(0,0.15), xlim=c(0.8,4.2), xaxt="n", cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3,  ylab=expression(paste("Average support 65%, ", delta[c1], "=0.4")), cex.lab=1.5, cex.axis=1.5, cex.main=1.5, main="", xlab="")
	segmentsL(0.98,mae.jo.0500.B04[1],0.98,mae.jo.0500.B04[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02, mae.ma.0500.B04[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02,mae.ma.0500.B04[1],1.02,mae.ma.0500.B04[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1, mae.di.0500.B04[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1,mae.di.0500.B04[1],1,mae.di.0500.B04[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 1
	points(0.98+q, mae.jo.0502.B04[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.0502.B04[1],0.98+q,mae.jo.0502.B04[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.0502.B04[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.0502.B04[1],1.02+q,mae.ma.0502.B04[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.0502.B04[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.0502.B04[1],1+q,mae.di.0502.B04[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 2
	points(0.98+q, mae.jo.0504.B04[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.0504.B04[1],0.98+q,mae.jo.0504.B04[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.0504.B04[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.0504.B04[1],1.02+q,mae.ma.0504.B04[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.0504.B04[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.0504.B04[1],1+q,mae.di.0504.B04[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 3
	points(0.98+q, mae.jo.0506.B04[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.0506.B04[1],0.98+q,mae.jo.0506.B04[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.0506.B04[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.0506.B04[1],1.02+q,mae.ma.0506.B04[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.0506.B04[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.0506.B04[1],1+q,mae.di.0506.B04[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	# lines
	segments(1,mae.jo.0500.B04[2],2,mae.jo.0502.B04[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(2,mae.jo.0502.B04[2],3,mae.jo.0504.B04[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(3,mae.jo.0504.B04[2],4,mae.jo.0506.B04[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(1,mae.ma.0500.B04[2],2,mae.ma.0502.B04[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(2,mae.ma.0502.B04[2],3,mae.ma.0504.B04[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(3,mae.ma.0504.B04[2],4,mae.ma.0506.B04[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(1,mae.di.0500.B04[2],2,mae.di.0502.B04[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(2,mae.di.0502.B04[2],3,mae.di.0504.B04[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(3,mae.di.0504.B04[2],4,mae.di.0506.B04[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))




axis(1, at=c(1,2,3,4), labels=c(expression(paste(rho,"=0.0")), expression(paste(rho,"=0.2")), expression(paste(rho,"=0.4")), expression(paste(rho,"=0.6")) ), lwd.ticks=0, cex.axis=1.5)

# second plot
q <- 0
plot(0.98, mae.jo.1000.B04[2], ylim=c(0,0.15), xlim=c(0.8,4.2), xaxt="n", cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3, ylab="", cex.lab=1.5, cex.axis=1.5, cex.main=1.5, main="", xlab="")
	segmentsL(0.98,mae.jo.1000.B04[1],0.98,mae.jo.1000.B04[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02, mae.ma.1000.B04[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02,mae.ma.1000.B04[1],1.02,mae.ma.1000.B04[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1, mae.di.1000.B04[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1,mae.di.1000.B04[1],1,mae.di.1000.B04[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 1
	points(0.98+q, mae.jo.1002.B04[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.1002.B04[1],0.98+q,mae.jo.1002.B04[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.1002.B04[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.1002.B04[1],1.02+q,mae.ma.1002.B04[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.1002.B04[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.1002.B04[1],1+q,mae.di.1002.B04[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 2
	points(0.98+q, mae.jo.1004.B04[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.1004.B04[1],0.98+q,mae.jo.1004.B04[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.1004.B04[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.1004.B04[1],1.02+q,mae.ma.1004.B04[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.1004.B04[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.1004.B04[1],1+q,mae.di.1004.B04[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 3
	points(0.98+q, mae.jo.1006.B04[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.1006.B04[1],0.98+q,mae.jo.1006.B04[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.1006.B04[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.1006.B04[1],1.02+q,mae.ma.1006.B04[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.1006.B04[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.1006.B04[1],1+q,mae.di.1006.B04[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	# lines
	segments(1,mae.jo.1000.B04[2],2,mae.jo.1002.B04[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(2,mae.jo.1002.B04[2],3,mae.jo.1004.B04[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(3,mae.jo.1004.B04[2],4,mae.jo.1006.B04[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(1,mae.ma.1000.B04[2],2,mae.ma.1002.B04[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(2,mae.ma.1002.B04[2],3,mae.ma.1004.B04[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(3,mae.ma.1004.B04[2],4,mae.ma.1006.B04[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(1,mae.di.1000.B04[2],2,mae.di.1002.B04[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(2,mae.di.1002.B04[2],3,mae.di.1004.B04[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(3,mae.di.1004.B04[2],4,mae.di.1006.B04[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))




axis(1, at=c(1,2,3,4), labels=c(expression(paste(rho,"=0.0")), expression(paste(rho,"=0.2")), expression(paste(rho,"=0.4")), expression(paste(rho,"=0.6")) ), lwd.ticks=0, cex.axis=1.5)



# third plot
q <- 0
plot(0.98, mae.jo.2000.B04[2], ylim=c(0,0.15), xlim=c(0.8,4.2), xaxt="n", cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3, ylab="", cex.lab=1.5, cex.axis=1.5, cex.main=1.5, main="", xlab="")
	segmentsL(0.98,mae.jo.2000.B04[1],0.98,mae.jo.2000.B04[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02, mae.ma.2000.B04[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02,mae.ma.2000.B04[1],1.02,mae.ma.2000.B04[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1, mae.di.2000.B04[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1,mae.di.2000.B04[1],1,mae.di.2000.B04[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 1
	points(0.98+q, mae.jo.2002.B04[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.2002.B04[1],0.98+q,mae.jo.2002.B04[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.2002.B04[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.2002.B04[1],1.02+q,mae.ma.2002.B04[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.2002.B04[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.2002.B04[1],1+q,mae.di.2002.B04[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 2
	points(0.98+q, mae.jo.2004.B04[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.2004.B04[1],0.98+q,mae.jo.2004.B04[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.2004.B04[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.2004.B04[1],1.02+q,mae.ma.2004.B04[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.2004.B04[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.2004.B04[1],1+q,mae.di.2004.B04[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 3
	points(0.98+q, mae.jo.2006.B04[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.2006.B04[1],0.98+q,mae.jo.2006.B04[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.2006.B04[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.2006.B04[1],1.02+q,mae.ma.2006.B04[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.2006.B04[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.2006.B04[1],1+q,mae.di.2006.B04[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	# lines
	segments(1,mae.jo.2000.B04[2],2,mae.jo.2002.B04[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(2,mae.jo.2002.B04[2],3,mae.jo.2004.B04[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(3,mae.jo.2004.B04[2],4,mae.jo.2006.B04[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(1,mae.ma.2000.B04[2],2,mae.ma.2002.B04[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(2,mae.ma.2002.B04[2],3,mae.ma.2004.B04[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(3,mae.ma.2004.B04[2],4,mae.ma.2006.B04[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(1,mae.di.2000.B04[2],2,mae.di.2002.B04[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(2,mae.di.2002.B04[2],3,mae.di.2004.B04[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(3,mae.di.2004.B04[2],4,mae.di.2006.B04[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))




axis(1, at=c(1,2,3,4), labels=c(expression(paste(rho,"=0.0")), expression(paste(rho,"=0.2")), expression(paste(rho,"=0.4")), expression(paste(rho,"=0.6")) ), lwd.ticks=0, cex.axis=1.5)








q <- 0
plot(0.98, mae.jo.0500.B085[2], ylim=c(0,0.15), xlim=c(0.8,4.2), xaxt="n", cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3,  ylab=expression(paste("Average support 80%, ", delta[c1], "=0.85")), cex.lab=1.5, cex.axis=1.5, cex.main=1.5, main="", xlab="")
	segmentsL(0.98,mae.jo.0500.B085[1],0.98,mae.jo.0500.B085[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02, mae.ma.0500.B085[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02,mae.ma.0500.B085[1],1.02,mae.ma.0500.B085[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1, mae.di.0500.B085[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1,mae.di.0500.B085[1],1,mae.di.0500.B085[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 1
	points(0.98+q, mae.jo.0502.B085[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.0502.B085[1],0.98+q,mae.jo.0502.B085[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.0502.B085[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.0502.B085[1],1.02+q,mae.ma.0502.B085[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.0502.B085[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.0502.B085[1],1+q,mae.di.0502.B085[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 2
	points(0.98+q, mae.jo.0504.B085[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.0504.B085[1],0.98+q,mae.jo.0504.B085[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.0504.B085[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.0504.B085[1],1.02+q,mae.ma.0504.B085[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.0504.B085[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.0504.B085[1],1+q,mae.di.0504.B085[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 3
	points(0.98+q, mae.jo.0506.B085[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.0506.B085[1],0.98+q,mae.jo.0506.B085[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.0506.B085[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.0506.B085[1],1.02+q,mae.ma.0506.B085[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.0506.B085[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.0506.B085[1],1+q,mae.di.0506.B085[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	# lines
	segments(1,mae.jo.0500.B085[2],2,mae.jo.0502.B085[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(2,mae.jo.0502.B085[2],3,mae.jo.0504.B085[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(3,mae.jo.0504.B085[2],4,mae.jo.0506.B085[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(1,mae.ma.0500.B085[2],2,mae.ma.0502.B085[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(2,mae.ma.0502.B085[2],3,mae.ma.0504.B085[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(3,mae.ma.0504.B085[2],4,mae.ma.0506.B085[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(1,mae.di.0500.B085[2],2,mae.di.0502.B085[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(2,mae.di.0502.B085[2],3,mae.di.0504.B085[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(3,mae.di.0504.B085[2],4,mae.di.0506.B085[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))




axis(1, at=c(1,2,3,4), labels=c(expression(paste(rho,"=0.0")), expression(paste(rho,"=0.2")), expression(paste(rho,"=0.4")), expression(paste(rho,"=0.6")) ), lwd.ticks=0, cex.axis=1.5)

# second plot
q <- 0
plot(0.98, mae.jo.1000.B085[2], ylim=c(0,0.15), xlim=c(0.8,4.2), xaxt="n", cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3, ylab="", cex.lab=1.5, cex.axis=1.5, cex.main=1.5, main="", xlab="")
	segmentsL(0.98,mae.jo.1000.B085[1],0.98,mae.jo.1000.B085[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02, mae.ma.1000.B085[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02,mae.ma.1000.B085[1],1.02,mae.ma.1000.B085[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1, mae.di.1000.B085[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1,mae.di.1000.B085[1],1,mae.di.1000.B085[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 1
	points(0.98+q, mae.jo.1002.B085[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.1002.B085[1],0.98+q,mae.jo.1002.B085[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.1002.B085[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.1002.B085[1],1.02+q,mae.ma.1002.B085[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.1002.B085[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.1002.B085[1],1+q,mae.di.1002.B085[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 2
	points(0.98+q, mae.jo.1004.B085[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.1004.B085[1],0.98+q,mae.jo.1004.B085[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.1004.B085[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.1004.B085[1],1.02+q,mae.ma.1004.B085[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.1004.B085[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.1004.B085[1],1+q,mae.di.1004.B085[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 3
	points(0.98+q, mae.jo.1006.B085[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.1006.B085[1],0.98+q,mae.jo.1006.B085[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.1006.B085[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.1006.B085[1],1.02+q,mae.ma.1006.B085[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.1006.B085[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.1006.B085[1],1+q,mae.di.1006.B085[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	# lines
	segments(1,mae.jo.1000.B085[2],2,mae.jo.1002.B085[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(2,mae.jo.1002.B085[2],3,mae.jo.1004.B085[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(3,mae.jo.1004.B085[2],4,mae.jo.1006.B085[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(1,mae.ma.1000.B085[2],2,mae.ma.1002.B085[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(2,mae.ma.1002.B085[2],3,mae.ma.1004.B085[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(3,mae.ma.1004.B085[2],4,mae.ma.1006.B085[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(1,mae.di.1000.B085[2],2,mae.di.1002.B085[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(2,mae.di.1002.B085[2],3,mae.di.1004.B085[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(3,mae.di.1004.B085[2],4,mae.di.1006.B085[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))




axis(1, at=c(1,2,3,4), labels=c(expression(paste(rho,"=0.0")), expression(paste(rho,"=0.2")), expression(paste(rho,"=0.4")), expression(paste(rho,"=0.6")) ), lwd.ticks=0, cex.axis=1.5)



# third plot
q <- 0
plot(0.98, mae.jo.2000.B085[2], ylim=c(0,0.15), xlim=c(0.8,4.2), xaxt="n", cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3, ylab="", cex.lab=1.5, cex.axis=1.5, cex.main=1.5, main="", xlab="")
	segmentsL(0.98,mae.jo.2000.B085[1],0.98,mae.jo.2000.B085[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02, mae.ma.2000.B085[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02,mae.ma.2000.B085[1],1.02,mae.ma.2000.B085[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1, mae.di.2000.B085[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1,mae.di.2000.B085[1],1,mae.di.2000.B085[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 1
	points(0.98+q, mae.jo.2002.B085[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.2002.B085[1],0.98+q,mae.jo.2002.B085[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.2002.B085[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.2002.B085[1],1.02+q,mae.ma.2002.B085[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.2002.B085[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.2002.B085[1],1+q,mae.di.2002.B085[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 2
	points(0.98+q, mae.jo.2004.B085[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.2004.B085[1],0.98+q,mae.jo.2004.B085[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.2004.B085[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.2004.B085[1],1.02+q,mae.ma.2004.B085[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.2004.B085[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.2004.B085[1],1+q,mae.di.2004.B085[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	q <- 3
	points(0.98+q, mae.jo.2006.B085[2], cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	segmentsL(0.98+q,mae.jo.2006.B085[1],0.98+q,mae.jo.2006.B085[3], col= rgb(255,0,0,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1.02+q, mae.ma.2006.B085[2], cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	segmentsL(1.02+q,mae.ma.2006.B085[1],1.02+q,mae.ma.2006.B085[3], col= rgb(0,0,255,100, maxColorValue=255), widL=0.08, lwd=2, lty=1)
	points(1+q, mae.di.2006.B085[2], cex=2.5, col= rgb(190,190,190,180, maxColorValue=255), pch=15)
	segmentsL(1+q,mae.di.2006.B085[1],1+q,mae.di.2006.B085[3], col= rgb(190,190,190,180, maxColorValue=255), widL=0.08, lwd=2, lty=1)

	# lines
	segments(1,mae.jo.2000.B085[2],2,mae.jo.2002.B085[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(2,mae.jo.2002.B085[2],3,mae.jo.2004.B085[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(3,mae.jo.2004.B085[2],4,mae.jo.2006.B085[2], lwd=1, col=rgb(255,0,0,100, maxColorValue=255))
	segments(1,mae.ma.2000.B085[2],2,mae.ma.2002.B085[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(2,mae.ma.2002.B085[2],3,mae.ma.2004.B085[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(3,mae.ma.2004.B085[2],4,mae.ma.2006.B085[2], lwd=1, col=rgb(0,0,255,100, maxColorValue=255))
	segments(1,mae.di.2000.B085[2],2,mae.di.2002.B085[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(2,mae.di.2002.B085[2],3,mae.di.2004.B085[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))
	segments(3,mae.di.2004.B085[2],4,mae.di.2006.B085[2], lwd=1, col=rgb(190,190,190,100, maxColorValue=255))




axis(1, at=c(1,2,3,4), labels=c(expression(paste(rho,"=0.0")), expression(paste(rho,"=0.2")), expression(paste(rho,"=0.4")), expression(paste(rho,"=0.6")) ), lwd.ticks=0, cex.axis=1.5)

#plot(1,1,col="white", bty="n", xaxt="none", yaxt="none", ylab="", xlab="", ylim=c(0,0.18), xlim=c(0.8,4.2))

plot(1,1,col="white", bty="n", xaxt="none", yaxt="none", ylab="", xlab="", ylim=c(0,0.18), xlim=c(0.8,4.2))
	points(1, 0.13, cex=2.5, col= "black", pch=19, lwd=3)
	segmentsL(1, 0.15,1,0.11, col= "black", widL=0.08, lwd=2, lty=1)
	text(1.05,0.13, "Mean of 100 Simulations", pos=4, cex=1.3)
	text(1.05,0.15, "0.975 Quantile", pos=4, cex=1.3)
	text(1.05,0.11, "0.025 Quantile", pos=4, cex=1.3)

plot(1,1,col="white", bty="n", xaxt="none", yaxt="none", ylab="", xlab="", ylim=c(0,0.18), xlim=c(0.8,4.2))


plot(1,1,col="white", bty="n", xaxt="none", yaxt="none", ylab="", xlab="", ylim=c(0,0.18), xlim=c(0.8,4.2))
	points(1, .15, cex=2.5, col= rgb(255,0,0,100, maxColorValue=255), pch=2, lwd=3)
	text(1.05, 0.15, "MrP Results", pos=4, cex=1.3)
	points(1, 0.13, cex=2.5, col= rgb(0,0,255,100, maxColorValue=255), pch=19)
	text(1.05, 0.13, "MrsP Results", pos=4, cex=1.3)
	points(1, 0.11, cex=2.5, col= rgb(190,190,190,100, maxColorValue=255), pch=15, lwd=3)
	text(1.05, 0.11, "Disaggregation Results", pos=4, cex=1.3)



dev.off()


